% Created 2021, Mike Siedlik
% Free for research and non-commercial use
% Please cite Siedlik and Issadore, 2021

% This MATLAB code was created with the following version of MATLAB:
% MATLAB Version: 9.10.0.1739362 (R2021a) Update 5
% This code uses the following MATLAB toolboxes:
% Image Processing Toolbox

% ---------------------------------------------------------
% ---------------------------------------------------------
% SPECIFY DROPLET DATA
% File name of image stack to be analyzed. This dataset is analogous to
% that created in A_ReconstructVideo, with the image centered on the moving
% droplet.
pic_name = 'C_Example.tif';
% File name of image stack of masks (roughly) corresponding to droplets
% listed above. These contain "edge regions" of the droplet that, when
% multiplied with the droplet image, correspond to dark pixel regions but
% don't correspond to washing. More representative mask images (that don't
% contain the dark edges) could be obtained by algorithms such as
% intensity-based thresholding.
mask_name = 'C_ExampleMasks.tif';

% Obtain useful dimensions of image data
INFO = imfinfo(pic_name);
im_dim = [INFO(1).Height,INFO(1).Width,numel(INFO)];


% Load reference droplet data. Here, load a data structure comparable to
% that obtained with B_CreateControlImages
load('C_cntl_data.mat');

% For each reference frame, crop out the regions close to the wall to avoid
% optical aberations that occur due to edge effects
for k = 1:5
    cntl_data(k).meanGreyValues = cntl_data(k).meanGreyValues(:,89:189);
    cntl_data(k).meanGreyValues([1:10,46:end],:) = 0;
end

% Rearrange the control data so that is in a form easily amenable to
% pixel-by-pixel linear regression
% Pixel grey values
pxVal = single(reshape([cntl_data.meanGreyValues],50,101,5));
% Known dilution
cVal = reshape([cntl_data.dilution],1,5);

% Matrix to store coefficients for linear regression model
M = zeros(size(pxVal,1:2));
% Matrix to store constants for linear regression model. Row/column indices
% correspond to 
B = zeros(size(pxVal,1:2));

% Create a separate linear model for each pixel
for r = 1:im_dim(1)
    for c = 1:im_dim(2)
        
        % If the pixel will not be considered, because it is too close to
        % an edge, don't create a model.
        if pxVal(r,c) == 0
            M(r,c) = nan;
            B(r,c) = nan;
        % Otherwise, do linear regression
        else
            C = [squeeze(pxVal(r,c,:)),ones(5,1)]\(cVal');
            M(r,c) = C(1);
            B(r,c) = C(2);
        end
    end
end

% Obtain current directory
DIR = cd;

% Create folder to save images
saveFolder = 'qMaps';
mkdir(saveFolder);

% Analyze each frame
for k = 1:im_dim(3)
    % 1) Read k-th frame of droplet stack
    dropPic = double(imread(pic_name,k));
    
    % 2) Read corresponding mask image
    maskPic = imread(mask_name,k);
    
    % 3) Take the product of 1) and 2) so that only internal droplet
    % regions will be operated on
    d2 = (dropPic.*(maskPic>0));
    d2(d2==0) = nan;
    
    % Perform pixel by pixel linear regression
    qMap = d2.*M + B;
    
    % Create plot
    f = figure('Visible','off');
    colormap(jet);
    surface(qMap,'LineStyle','none');
    
    % Plot formatting
    axis image
    cbar = colorbar;
    caxis([0,1]);
    ylim([1,50]);
    xlim([1,101]);
    set(gca,'CLim',[0,0.9]);
    cbar.Ticks = 0:0.2:1;
    cbar.TickLabels = num2cell(0:0.2:1);
    cbar.Limits = [0,1];
    set(gca,'YDir','reverse');
    set(gca,'Units','pixels');
    set(gca,'YTick','');
    set(gca,'XTick','');
    set(gca,'Color','k');
    set(gca,'YColor','k');
    set(gca,'XColor','k');
    cbar.Label.String = 'C/C_0';
    cbar.Label.FontSize = 12;
    set(gca,'Box','on');
    
    % Print image to save folder
    cd(strcat(DIR,'\',saveFolder));
    print(strcat(num2str(k),'.tif'),'-dtiff');
    cd(DIR);
    
    close all
    disp(k)
end



